
## Clear environment and return memory to the operating system
rm(list = ls())
gc()

## Load (install) required packages
suppressMessages({
  library(tidyverse)
  library(commarobust)
  library(xtable)
  library(ggthemes)
  library(cowplot)
  library(MBESS)
  library(ri2)
  library(estimatr)
  library(grf)
  library(scales)
  })

## Helper functions
add_parens <- function(x, digits = 2) {
  x <- as.numeric(x)
  return(paste0("(", sprintf(paste0("%.", digits, "f"), x), ")"))
}

format_num <- function(x, digits = 2) {
  x <- as.numeric(x)
  return(paste0(sprintf(paste0("%.", digits, "f"), x)))
}

make_entry <- function(est, se, p, digits = 2) {
  entry <- paste0(format_num(est, digits = digits), " ", 
                  add_parens(se, digits = digits))
  entry[p < 0.05] <- paste0(entry[p < 0.05], "*")
  return(entry)
}

table_entry <- function(est, se, digits = 2) {
  entry <- paste0(format_num(est, digits = digits), " ", 
                  add_parens(se, digits = digits))
  return(entry)
}

standardize <- function(x){
  (x-min(x, na.rm = TRUE))/(max(x, na.rm = TRUE) - min(x, na.rm = TRUE))
}

glass_delta <- function(Y, Z, name = "Placebo"){
  Y/sd(Y[Z == name], na.rm = TRUE)
}

## Theme for pretty plots
devtools::source_url("https://kylepeyton.github.io/assets/theme_kyle.R")

## Make some adjustments ...
theme_kyle2 <- function () { 
  theme_kyle(font_size = 22) + #%+replace%
    theme(legend.position = "bottom",
          legend.justification = "center",
          text = element_text(family = "serif", size = 18),
          plot.title = element_text(margin = margin(b = 0, l = -10), 
                                    size = 24, face = "bold", family = "serif"), 
          legend.key.size = unit(3, "line"),
          legend.margin = margin(0,0,0,0),
          legend.box.margin = margin(-30, -10, -10, -10),
          strip.text.x = element_text(size = 22, family = "serif"),
          axis.title.y = element_text(size = 19, family = "serif"),
          legend.title = element_blank(),
          axis.text.x = element_text(size = 19, family = "serif"),
          strip.background = element_rect(fill = "grey90"),
          axis.line.y = element_line(size = I(4/12), colour = "black"),
          axis.ticks = element_line(colour = "black"),
          axis.ticks.y = element_line(size = I(4/12), colour = "black"))
}

## Set paths for where data is stored, where figures/tables should be exported
datapath <- "/XXX/Data/"
plotpath <- "/XXX/Figures/"
tabpath <- "/XXX/Tables/"

## Read in dataset 
combined_trust <- read.csv(paste0(datapath, "combined_trust.csv"))

## Construct metric instrumental variable (Z is vector of random assignments)
combined_trust <- 
  combined_trust %>%
  mutate(Z = factor(Z, levels = c("Corrupt", "Placebo", "Honest")),
         Z1 = NA, 
         Z1 = ifelse(Z == "Honest", 1, Z1),
         Z1 = ifelse(Z == "Placebo", 0.5, Z1),
         Z1 = ifelse(Z == "Corrupt", 0, Z1))

## Studies ordered by date: Experiment 1-2, Experiment 4-5, Experiment 3
study <- c("MTurk_Pols14", "GenPop_Pols14", "MTurk_NFL14", "MTurk_NFL15",
           "MTurk_Pols17")

## Create dataframe for estimation and dataframe for placebo experiments. 
## Standardize DVs and trust measures
est_df <- 
  combined_trust %>% 
  ungroup() %>% 
  filter(experiment %in% study[c(1, 2, 5)]) %>% 
  mutate(D = glass_delta(Y = fmptrust_index, Z = Z),
         D1 = glass_delta(Y = trustgov_scale, Z = Z),
         Y = glass_delta(Y = fedspend_index, Z = Z))

placebo_df <- 
  combined_trust %>% 
  ungroup() %>% 
  filter(experiment %in% study[c(3, 4)]) %>% 
  mutate(D = glass_delta(Y = fmptrust_index, Z = Z),
         D1 = glass_delta(Y = trustgov_scale, Z = Z),
         Y = glass_delta(Y = fedspend_index, Z = Z))

#####--------------------------Figure 2------------------------------------#####
gg_df <- 
  combined_trust %>% 
  group_by(study_type) %>% 
  mutate(D = glass_delta(Y = fmptrust_index, Z = Z),
         Y = glass_delta(Y = fedspend_index, Z = Z)) %>% 
  ungroup() %>% 
  mutate(Z = factor(Z, levels = c("Corrupt", "Placebo", "Honest"),
                    labels = c("Corrupt", "Control", "Honest")),
         experiment = factor(experiment, 
                             levels = c("MTurk_Pols14", "GenPop_Pols14",
                                        "MTurk_Pols17", "MTurk_NFL14",
                                        "MTurk_NFL15"),
                             labels = c("MTurk June '14 \n (N = 643) ",
                                        "National Sept. '14 \n (N = 1,324)",
                                        "MTurk Mar. '17 \n (N = 1,870)",
                                        "MTurk Dec. '14 \n (N = 584)",
                                        "MTurk July '15 \n (N = 585)")),
         study_group = ifelse(study_type == "NFL", "Placebo Experiments",
                              "Corruption Experiments")) 

avg_df <- 
  gg_df %>% 
  group_by(Z, experiment) %>% 
  do(tidy(lm_robust(fmptrust_scale ~ 1, weights = wts, data = .))) %>% 
  select(-term) %>% 
  mutate(measure = "Trust in Government",
         study_group = ifelse(experiment %in% 
                                c("MTurk June '14 \n (N = 643) ",
                                  "National Sept. '14 \n (N = 1,324)",
                                  "MTurk Mar. '17 \n (N = 1,870)"),
                              "Corruption Experiments", 
                              "Placebo Experiments")) %>%
  rbind(gg_df  %>%
          group_by(Z, experiment) %>% 
          do(tidy(lm_robust(fedspend_scale ~ 1, weights = wts, data = .))) %>% 
          select(-term) %>%  
          mutate(measure = "Support for Redistribution",
                 study_group = ifelse(experiment %in% 
                                        c("MTurk June '14 \n (N = 643) ",
                                          "National Sept. '14 \n (N = 1,324)",
                                          "MTurk Mar. '17 \n (N = 1,870)"),
                                      "Corruption Experiments", 
                                      "Placebo Experiments"))) %>% 
  mutate(measure = factor(measure, levels = c("Trust in Government",
                                              "Support for Redistribution")),
         group_lab = factor(
           paste(measure, study_group, sep = "\n"),
           levels = c("Trust in Government\nCorruption Experiments",
                      "Trust in Government\nPlacebo Experiments",
                      "Support for Redistribution\nCorruption Experiments",
                      "Support for Redistribution\nPlacebo Experiments")))

pooled_df <- 
  gg_df %>% 
  group_by(Z, study_group) %>% 
  do(tidy(lm_robust(D ~ 1, weights = wts, data = .))) %>% 
  filter(term == "(Intercept)") %>% 
  select(-term) %>% 
  mutate(measure = "Trust in Government") %>% 
  rbind(gg_df %>% 
          group_by(Z, study_group) %>% 
          do(tidy(lm_robust(Y ~ 1, weights = wts, data = .))) %>% 
          filter(term == "(Intercept)") %>% 
          select(-term) %>% 
          mutate(measure = "Support for Redistribution")) %>% 
  mutate(measure = factor(measure, levels = c("Trust in Government",
                                              "Support for Redistribution")),
         group_lab = factor(
           paste(measure, study_group, sep = "\n"),
           levels = c("Trust in Government\nCorruption Experiments",
                      "Trust in Government\nPlacebo Experiments",
                      "Support for Redistribution\nCorruption Experiments",
                      "Support for Redistribution\nPlacebo Experiments")))

## Plot of First Stage Effects
first_stage_means_plot <- 
  pooled_df %>% filter(measure == "Trust in Government") %>% 
  ggplot(., aes(y = estimate, x = Z)) +
  geom_crossbar(aes(x = Z, y = estimate, ymin = conf.low, ymax = conf.high), 
                fatten = 1,
                color = "black", fill = "gray",
                alpha = 0.5, inherit.aes = FALSE) +
  facet_wrap(~ study_group) + 
  geom_linerange(data = avg_df %>% filter(measure == "Trust in Government"), 
                 aes(x = Z, color = experiment, ymax = conf.high, 
                     ymin = conf.low), inherit.aes = FALSE, show.legend = FALSE,
                 position = position_dodge(width = 0.9)) +
  geom_point(data = avg_df %>% filter(measure == "Trust in Government"), 
             aes(y = estimate, x = Z, color = experiment, shape = experiment,
                 size = experiment),
             position = position_dodge(width = 0.9),
             stroke = 1, inherit.aes = FALSE, show.legend = TRUE)  +
  scale_color_grey("", start = 0, end = 0) +
  scale_size_manual("Sample:", values = c(2.5, 2.5, 3, 2.5, 2.5)) + 
  scale_shape_manual("Sample:", values = c(15, 20, 18, 0, 1)) + 
  xlab("") + 
  ylab("Trust in Government\n (standardized scale)")  +
  theme_kyle2() + 
  guides(color = FALSE, size = FALSE,
         shape = guide_legend(order = 1, override.aes = list(size = 3),
                              nrow = 2, byrow = TRUE)) +
  scale_y_continuous(limits = c(1.7, 3.3), breaks = seq(1, 3.2, 0.2),
                     expand = c(0.01, 0)) +
  theme(
    axis.line.x = element_line(size = 0.1, colour = "grey70"),
    axis.line.y = element_line(size = 0.1, colour = "grey70"),
    legend.position = "bottom",
    legend.key.size = unit(1, "line"),
    legend.title = element_text(size = 10, family = "serif"),
    legend.text = element_text(size = 10, family = "serif"),
    legend.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "cm"),
    plot.margin = unit(c(t = 1, r = 0.5, b = 1, l = 0.5), "lines"),
    plot.title = element_text(hjust = 0),
    strip.text.x = element_text(size = 14, family = "serif"),
    axis.title = element_text(size = 14, family = "serif"),
    axis.title.y = element_text(size = 14, family = "serif"),
    axis.text.x = element_text(size = 14, family = "serif", face = "italic",
                               margin = margin(t = 0)),
    axis.text.y = element_text(size = 14, family = "serif"),
    axis.ticks.y = element_line(colour = "grey70"),
    axis.ticks.x = element_blank(),
    panel.grid.minor = element_blank(),
    panel.grid.major.y  = element_blank(),
    panel.border = element_rect(fill = NA, colour = "grey70")
  ) 

ggsave(first_stage_means_plot, dpi = 1000,
       filename = paste0(plotpath, "first_stage_means_plot.pdf"),
       height = 4, width = 6, device = "pdf")
setwd(plotpath)
system("open first_stage_means_plot.pdf")

#####--------------------------Figure 3, Table S21-------------------------#####

# ITTd:
itt_d_fit1 <- 
  lm_robust(D ~ Z1 + dem_pid3 + dem_age + dem_female + dem_college + dem_work +
              dem_income + dem_white + dem_ideo5 + experiment, weights = wts,
            data = est_df) 

itt_d_fit2 <- 
  lm_robust(D ~ Z1 + dem_pid3 + dem_age + dem_female + dem_college + dem_work + 
              dem_income + dem_white + dem_ideo5 + experiment, weights = wts,
            data = placebo_df) 
# ITTy:
itt_y_fit1 <- 
  lm_robust(Y ~ Z1 + dem_pid3 + dem_age + dem_female + dem_college + dem_work + 
              dem_income + dem_white + dem_ideo5 + experiment, weights = wts,
            data = est_df) 

itt_y_fit2 <- 
  lm_robust(Y ~ Z1 + dem_pid3 + dem_age + dem_female + dem_college + dem_work + 
              dem_income + dem_white + dem_ideo5 + experiment, weights = wts,
            data = placebo_df) 

itt_d_est <- 
  itt_d_fit1 %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c('(Intercept)', "Z1")) %>% 
  mutate(S = "Political Corruption") %>% 
  bind_rows(
    itt_d_fit2 %>% 
      tidy() %>% 
      mutate_if(is.numeric, round, 3) %>% 
      filter(term %in% c('(Intercept)', "Z1")) %>% 
      mutate(S = "Non-Political Corruption")
  ) %>% 
  mutate(estimand = "Trust in Government (First Stage)")


itt_y_est <-
  itt_y_fit1 %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c('(Intercept)', "Z1")) %>% 
  mutate(S = "Political Corruption") %>% 
  bind_rows(
    itt_y_fit2 %>% 
      tidy() %>% 
      mutate_if(is.numeric, round, 3) %>% 
      filter(term %in% c('(Intercept)', "Z1")) %>% 
      mutate(S = "Non-Political Corruption")
  ) %>% 
  mutate(estimand = "Support for Redistribution (Reduced Form)")


itt_est_df <- bind_rows(itt_d_est, itt_y_est) %>% filter(term == "Z1")


itt_plot_df <- 
  itt_est_df %>% 
  mutate(S = ifelse(S == "Political Corruption", 
                    "\nCorruption\nExperiments\n (N = 3,837)", 
                    "Placebo\nExperiments\n (N = 1,169)"),
         S = factor(S, levels = rev(c("\nCorruption\nExperiments\n (N = 3,837)",
                                      "Placebo\nExperiments\n (N = 1,169)"))),
         estimand = factor(estimand, 
                           levels = c("Trust in Government (First Stage)",
                                      "Support for Redistribution (Reduced Form)"),
                           labels = c("Trust in Government \n (First Stage Effect)", 
                                      "Support for Redistribution \n (Reduced Form Effect)"))) 


itt_plot_equiv <- 
  itt_plot_df %>% 
  ggplot(., aes(x = estimate, y = S, color = S, group = S, shape = S)) + 
  geom_vline(xintercept = 0, size = 0.5, lty = 1, colour = "grey70") + 
  geom_vline(xintercept = 0.20, size = 0.5, lty = 2, colour = "grey70") + 
  geom_vline(xintercept = -0.20, size = 0.5, lty = 2, colour = "grey70") + 
  geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), height = 0,
                 size = 0.65, color = "black", alpha = 0.8) + 
  geom_errorbarh(aes(xmin = estimate - 1.65*std.error, 
                     xmax = estimate + 1.65*std.error), 
                 size = 1.25, color = "black", height = 0) +
  geom_point(size = 3.5, col = "black", fill = "black", show.legend = FALSE) +
  scale_shape_manual(" ", values = c(15, 16)) +
  ylab("") +   xlab("Covariate-Adjusted OLS Estimate (in Standard Units)") +
  facet_wrap(~ estimand) +
  scale_x_continuous(limits = c(-0.35, 0.75), breaks = seq(-0.20, 0.60, 0.20), 
                     labels = c("-0.2", "0.0", "0.2", "0.4", "0.6"),
                     expand = c(0, 0)) +
  scale_y_discrete(drop = TRUE) +
  theme_kyle2() + 
  theme(
    axis.line.x = element_line(size = 0.1, colour = "grey70"),
    axis.line.y = element_line(size = 0.1, colour = "grey70"),
    legend.position = "none",
    legend.key.size = unit(1, "line"),
    legend.title = element_text(size = 10, family = "serif"),
    legend.text = element_text(size = 10, family = "serif"),
    legend.margin = margin(t = 0, r = 0, b = 0, l = 0, unit = "cm"),
    plot.margin = unit(c(t = 0.2, r = 0.2, b = 1, l = 0.2), "lines"),
    strip.text.x = element_text(size = 14, family = "serif"),
    axis.title = element_text(size = 14, family = "serif"),
    axis.text.x = element_text(size = 12, family = "serif"),
    axis.text.y = element_text(size = 14, family = "serif"),
    axis.title.y = element_blank(),
    axis.ticks.x = element_line(colour = "grey70"),
    axis.ticks.y = element_blank(),
    panel.grid.minor = element_blank(),
    panel.grid.major.y  = element_blank(),
    panel.border = element_rect(fill = NA, colour = "grey70")
  )

ggsave(itt_plot_equiv,
       filename = paste0(plotpath, "itt_plot_equiv.pdf"),
       height = 4, width = 7, device = "pdf", dpi = 1000)
setwd(plotpath)
system("open itt_plot_equiv.pdf")

# Estimate LATE of D on Y in politics experiments
late_fit1 <- 
  iv_robust(Y ~ D + dem_pid3 + dem_age + dem_female + dem_college + dem_income + 
              dem_work + dem_white + dem_ideo5 + experiment | Z1 + dem_pid3 + 
              dem_age + dem_female + dem_college + dem_income + dem_work + 
              dem_white + dem_ideo5 + experiment, weights = wts,
            data = est_df) 

late_est_df <- 
  late_fit1 %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c('(Intercept)', "D")) %>% 
  mutate(estimand = "LATE", S = "Political Corruption") 

## Print estimates/p-values for reporting in manuscript
itt_plot_df %>% 
  mutate(ci_low_90 = estimate - 1.65*std.error,
         ci_hi_90 = estimate + 1.65*std.error) %>% 
  select(S, estimand, estimate, std.error, statistic, p.value, conf.low, 
         conf.high, ci_low_90, ci_hi_90) %>% 
  mutate_if(is.numeric, round, 2)

late_est_df %>% 
  filter(term == "D") %>% 
  mutate(ci_low_90 = estimate - 1.65*std.error,
         ci_hi_90 = estimate + 1.65*std.error) %>% 
  select(S, estimand, estimate, std.error, statistic, p.value, conf.low, 
         conf.high, ci_low_90, ci_hi_90) %>% 
  mutate_if(is.numeric, round, 2)

## Group means by partisanship in Placebo condition
combined_trust %>% 
  filter(Z == "Placebo", dem_pid3 != "Independent") %>% 
  group_by(dem_pid3) %>% 
  do(tidy(lm_robust(fedspend_scale ~ 1, weights = wts, data = .))) %>% 
  select(-term)

# Export table for SI: 
pooled_regtable <-
  bind_rows(itt_est_df, late_est_df) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(entry = make_entry(est = estimate, se = std.error, p = p.value)) %>% 
  select(estimand, entry, S) %>%
  spread(estimand, entry) %>%
  arrange(desc(S)) %>% 
  select(S, `Trust in Government (First Stage)`, 
         `Support for Redistribution (Reduced Form)`, LATE)

pooled_regtable[2,4] <- "-"

xtable(pooled_regtable) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "pooled_regtable.tex"))


#####--------------------------Figure 4------------------------------------#####

covs <- c("dem_pid3", "dem_age", "dem_female", "dem_college", "dem_income", 
          "dem_white", "dem_ideo5", "dem_work")

# NB: GRF only accepts rows with complete cases on DVs and all covariates
index <- with(est_df, complete.cases(fmptrust_index, fedspend_index))
X <- est_df[index, c(covs, "experiment")]
X$study1 <- as.numeric(X$experiment == "MTurk_Pols14")
X$study2 <- as.numeric(X$experiment == "GenPop_Pols14")
X$study3 <- as.numeric(X$experiment == "MTurk_Pols17")
X$dem_republican <- as.numeric(X$dem_pid3 == "Republican")
X$dem_democrat <- as.numeric(X$dem_pid3 == "Democrat")
X$dem_independent <- as.numeric(X$dem_pid3 == "Independent")
X <- X[, -which(names(X) %in% c("dem_pid3", "experiment"))]
X <- as.matrix(X)

Y_Obs <- est_df$fedspend_index[index]
W <- est_df$fmptrust_index[index]
Z <- est_df$Z1[index]

# Step 1. Train a causal forest using the observed outcome (Y_Obs), the matrix
# of pre-treatment covariates (X), and the treatment assignment vector (W).
set.seed(123)
tau_forest <- instrumental_forest(Y = Y_Obs, X = X, W = W, Z = Z, 
                                  honesty = TRUE, num.trees = 4000, seed = 123)

# Step 2. Estimate the treatment effects for all units using
# out-of-bag prediction.
tau_hat_forest <- predict(tau_forest, estimate.variance = TRUE)

# Step 3. Compute estimated SEs treatment effects for predictions. 
sigma_hat <- sqrt(tau_hat_forest$variance.estimates)

# Step 4. Store the estimated treatment effects and CIs in a data frame.
grf_df <- 
  cbind(data.frame(
    tau_hat = tau_hat_forest$predictions,
    ci_upper = tau_hat_forest$predictions + 1.96 * sigma_hat,
    ci_lower = tau_hat_forest$predictions - 1.96 * sigma_hat),
    X
  )

# Order treatment effects by size and generate plot. 
grf_df$order <- with(grf_df, rank(tau_hat, ties.method = "first"))
grf_df$study <- NA
grf_df$study[grf_df$study1 == 1] <- "Study 1"
grf_df$study[grf_df$study2 == 1] <- "Study 2"
grf_df$study[grf_df$study3 == 1] <- "Study 3"

plot_df <-
  grf_df %>%
  group_by(study) %>%
  mutate(order = rank(tau_hat, ties.method = "first")) %>%
  ungroup()

grf_plot <- 
  grf_df %>% 
  ggplot(., aes(x = tau_hat, y = order, color = study)) + 
  geom_errorbarh(aes(xmin = ci_lower, xmax = ci_upper), col = "dimgrey",
                 alpha = 0.2, height = 0, size = 0.5) + 
  geom_point(size = 1, alpha = 0.8, pch = 20, col = "black", fill = "black") +
  geom_vline(xintercept = 0, size = 0.5, lty = 2) +
  labs(title = "",
       x = "Instrumental Forest Estimated Treatment Effect", 
       y = "") + 
  theme_kyle2() + 
  theme(
    axis.line.x = element_line(size = 0.5, colour = "black"),
    axis.line.y = element_blank(),
    legend.position = "none",
    plot.margin = unit(c(t = -0.5, r = 0.4, b = 1, l = 0.2), "lines"),
    axis.title = element_text(size = 14, family = "serif"),
    axis.title.x = element_text(size = 14, family = "serif"),
    axis.text = element_text(size = 14, family = "serif"),
    axis.text.x = element_text(size = 14, family = "serif"),
    axis.title.y = element_blank(),
    axis.ticks.x = element_line(colour = "black"),
    axis.ticks.y = element_blank(),
    axis.text.y = element_blank(),
    panel.border = element_rect(fill = NA, colour = "black", size = 0.5)) +
  scale_x_continuous(limits = c(-0.45, 0.45), breaks = seq(-0.4, 0.4, 0.20)) +
  scale_y_continuous(expand = c(0.01, 0.01))


ggsave(grf_plot,
       filename = paste0(plotpath, "grf_plot.pdf"),
       height = 5, width = 4, device = "pdf", dpi = 1000)
setwd(plotpath)
system("open grf_plot.pdf")

# How many CIs include zero? 
grf_df %>% 
  summarise(mean(ci_lower <= 0 & ci_upper >= 0))

# What proportion of the CIs that don't include zero are positive? 
these <- with(grf_df, which(ci_lower <= 0 & ci_upper >= 0))
grf_df[-these,] %>% 
  filter(tau_hat > 0) %>% 
  nrow()/nrow(grf_df)


#####--------------------------Figures S2-S6-------------------------------#####

## Function for randomization inference (RI) for covariate balance: regress 
## instrument on background covariates and extract f-statistic
balance_fun <- function(data){
  summary(lm(as.integer(Z) ~ factor(dem_partyid) + dem_age + dem_female + 
               factor(dem_edu) + scale(dem_income) + dem_ideo5 + dem_work + 
               factor(dem_race5), data = data))$f[1]
}

set.seed(123) # for reproducability, run everything in this section concurrently

## RI for Experiment 1
exp1_df <- combined_trust %>% filter(experiment == "MTurk_Pols14")
design_exp1 <- declare_ra(N = nrow(exp1_df), num_arms = 3, simple = TRUE,
                          conditions = c("Placebo", "Corrupt", "Honest"))

ri_exp1 <- 
  conduct_ri(
    test_function = balance_fun,
    declaration = design_exp1, 
    assignment = "Z", 
    sharp_hypothesis = 0,
    data = exp1_df, 
    sims = 5000
  )

## RI for Experiment 2
exp2_df <- combined_trust %>% filter(experiment == "GenPop_Pols14")
block_mat <- as.matrix(table(exp2_df$dem_pid3, exp2_df$Z))
design_exp2 <- declare_ra(N = nrow(exp2_df), num_arms = 3, 
                          blocks = exp2_df$dem_pid3,
                          block_m_each = block_mat,
                          conditions = c("Placebo", "Corrupt", "Honest"))

ri_exp2 <- 
  conduct_ri(
    test_function = balance_fun,
    declaration = design_exp2, 
    assignment = "Z", 
    sharp_hypothesis = 0,
    data = exp2_df, 
    IPW_weights = "wts",
    sims = 5000
  )

## RI for Experiment 3
exp3_df <- combined_trust %>% filter(experiment == "MTurk_Pols17")
design_exp3 <- declare_ra(N = nrow(exp3_df), num_arms = 3, simple = TRUE,
                          conditions = c("Placebo", "Corrupt", "Honest"))

ri_exp3 <- 
  conduct_ri(
    test_function = balance_fun,
    declaration = design_exp3, 
    assignment = "Z", 
    sharp_hypothesis = 0,
    data = exp3_df,
    sims = 5000
  )



## RI for Experiment 4
exp4_df <- combined_trust %>% filter(experiment == "MTurk_NFL14")
design_exp4 <- declare_ra(N = nrow(exp4_df), num_arms = 3, simple = TRUE,
                          conditions = c("Placebo", "Corrupt", "Honest"))


ri_exp4 <- 
  conduct_ri(
    test_function = balance_fun,
    declaration = design_exp4, 
    assignment = "Z", 
    sharp_hypothesis = 0,
    data = exp4_df,
    sims = 5000
  )

## RI for Experiment 5
exp5_df <- combined_trust %>% filter(experiment == "MTurk_NFL15")
design_exp5 <- declare_ra(N = nrow(exp5_df), num_arms = 3, simple = TRUE,
                          conditions = c("Placebo", "Corrupt", "Honest"))

ri_exp5 <- 
  conduct_ri(
    test_function = balance_fun,
    declaration = design_exp5, 
    assignment = "Z", 
    sharp_hypothesis = 0,
    data = exp5_df,
    sims = 5000
  )


## Summarise results and export figures
ri_out <- bind_rows(summary(ri_exp1), summary(ri_exp2),summary(ri_exp3),
                    summary(ri_exp4), summary(ri_exp5)) %>%
  mutate(experiment = paste("Experiment", 1:5))

ri_out %>% mutate_if(is.numeric, round, 2)

ri_exp1$sims_df$term <- "F-statistic"
p <- plot(ri_exp1) + ggtitle("")

ggsave(paste0(plotpath, "exp1_ri_plot.pdf"), p, 
       device = "pdf", width = 5, height = 5)
setwd(plotpath)
system("open exp1_ri_plot.pdf")

ri_exp2$sims_df$term <- "F-statistic"
p <- plot(ri_exp2) + ggtitle("")

ggsave(paste0(plotpath, "exp2_ri_plot.pdf"), p, 
       device = "pdf", width = 5, height = 5)
setwd(plotpath)
system("open exp2_ri_plot.pdf")

ri_exp3$sims_df$term <- "F-statistic"
p <- plot(ri_exp3) + ggtitle("")

ggsave(paste0(plotpath, "exp3_ri_plot.pdf"), p, 
       device = "pdf", width = 5, height = 5)
setwd(plotpath)
system("open exp3_ri_plot.pdf")

ri_exp4$sims_df$term <- "F-statistic"
p <- plot(ri_exp4) + ggtitle("")

ggsave(paste0(plotpath, "exp4_ri_plot.pdf"), p, 
       device = "pdf", width = 5, height = 5)
setwd(plotpath)
system("open exp4_ri_plot.pdf")

ri_exp5$sims_df$term <- "F-statistic"
p <- plot(ri_exp5) + ggtitle("")

ggsave(paste0(plotpath, "exp5_ri_plot.pdf"), p, 
       device = "pdf", width = 5, height = 5)
setwd(plotpath)
system("open exp5_ri_plot.pdf")


#####--------------------------Figures S7-S8-------------------------------#####
anes_trustgov_plot <- 
  combined_trust %>%
  filter(Z == "Placebo") %>%
  mutate(study = NA, 
         study = ifelse(experiment == "MTurk_Pols14", "Experiment 1", 
                        study),
         study = ifelse(experiment == "GenPop_Pols14", "Experiment 2", 
                        study),
         study = ifelse(experiment == "MTurk_NFL14", "Experiment 4", 
                        study),
         study = ifelse(experiment == "MTurk_NFL15", "Experiment 5", 
                        study),
         study = ifelse(experiment == "MTurk_Pols17", "Experiment 3", 
                        study),
         study = factor(study, levels = c("Experiment 1", 
                                          "Experiment 2", 
                                          "Experiment 3",
                                          "Experiment 4", 
                                          "Experiment 5"))) %>%
  ggplot(., aes(x = trustgov_trustgstd)) + 
  geom_bar(aes(y = (..count..)/tapply(..count..,..PANEL..,sum)[..PANEL..])) + 
  facet_wrap(~ study, scales = "free_x", ncol = 2) + 
  scale_y_continuous(labels = percent_format()) + 
  xlab("") + ylab("") +
  theme_kyle2()

ggsave(paste0(plotpath, "anes_control_raw.pdf"), anes_trustgov_plot, 
       device = "pdf", width = 7, height = 10)
setwd(plotpath)
system("open anes_control_raw.pdf")


fmp_trustgov_plot <- 
  combined_trust %>%
  filter(Z == "Placebo") %>%
  mutate(study = NA, 
         study = ifelse(experiment == "MTurk_Pols14", "Experiment 1", 
                        study),
         study = ifelse(experiment == "GenPop_Pols14", "Experiment 2", 
                        study),
         study = ifelse(experiment == "MTurk_NFL14", "Experiment 4", 
                        study),
         study = ifelse(experiment == "MTurk_NFL15", "Experiment 5", 
                        study),
         study = ifelse(experiment == "MTurk_Pols17", "Experiment 3", 
                        study),
         study = factor(study, levels = c("Experiment 1", 
                                          "Experiment 2", 
                                          "Experiment 3",
                                          "Experiment 4", 
                                          "Experiment 5"))) %>%
  filter(!is.na(fmptrust_index)) %>% 
  ggplot(., aes(x = fmptrust_index)) + 
  geom_bar(aes(y = (..count..)/tapply(..count..,..PANEL..,sum)[..PANEL..])) + 
  facet_wrap(~ study, ncol = 2) + 
  scale_y_continuous(labels = percent_format()) + 
  scale_x_continuous(breaks = seq(4, 24, by = 4)) + 
  xlab("") + ylab("") +
  theme_kyle2()


ggsave(paste0(plotpath, "fmp_control_raw.pdf"), fmp_trustgov_plot, 
       device = "pdf", width = 7, height = 10)
setwd(plotpath)
system("open fmp_control_raw.pdf")


#####--------------------------Figures S13-S14----------------------------#####
post_trend_plot <- 
  combined_trust %>% 
  filter(experiment == "MTurk_Pols17" & !is.na(post_direction)) %>%
  mutate(Z = factor(Z, 
                    levels = c("Corrupt", "Placebo", "Honest"), 
                    labels = c("Corrupt", "Control", "Honest"))) %>%
  ggplot(., aes(x = factor(post_direction))) + 
  geom_bar(aes(y = (..count..)/tapply(..count..,..PANEL..,sum)[..PANEL..])) + 
  scale_x_discrete(labels = c("Strongly \n decreased", "Somewhat \n decreased",
                              "Neither increased \n nor decreased", 
                              "Somewhat \n increased", "Strongly \n increased")) + 
  facet_wrap(~ Z) + 
  scale_y_continuous(labels = percent_format()) + 
  xlab("") + ylab("") +
  theme_kyle2() +
  theme(axis.text.x = element_text(angle = 40, hjust = 1, size = 10)) 

post_support_plot <- 
  combined_trust %>% 
  filter(experiment == "MTurk_Pols17" & !is.na(post_support)) %>%
  mutate(Z = factor(Z, 
                    levels = c("Corrupt", "Placebo", "Honest"), 
                    labels = c("Corrupt", "Control", "Honest"))) %>%
  ggplot(., aes(x = factor(post_support))) + 
  geom_bar(aes(y = (..count..)/tapply(..count..,..PANEL..,sum)[..PANEL..])) + 
  scale_x_discrete(labels = rev(c("Strongly \n supports", "Somewhat \n supports",
                                  "Neither supports \n nor undermines", 
                                  "Somewhat \n undermines", "Strongly \n undermines"))) + 
  facet_wrap(~ Z) + 
  scale_y_continuous(labels = percent_format()) + 
  xlab("") + ylab("") +
  theme_kyle2() + 
  theme(axis.text.x = element_text(angle = 40, hjust = 1, size = 10)) 

ggsave(paste0(plotpath, "post_trend_plot.png"), post_trend_plot, 
       device = "png", width = 12, height = 8)
setwd(plotpath)
system("open post_trend_plot.png")

ggsave(paste0(plotpath, "post_support_plot.png"), post_support_plot, 
       device = "png", width = 12, height = 8)
setwd(plotpath)
system("open post_support_plot.png")


#####--------------------------Tables S2-S13------------------------------#####

covs <- c("dem_female", "dem_age", "dem_college","dem_work", "dem_asian", 
          "dem_black", "dem_hispanic", "dem_white", "dem_democrat", 
          "dem_republican")


for(i in 1:length(study)){
  bal <- with(combined_trust %>% filter(experiment == study[i]), table(Z))
  bal <- c(bal, sum(bal))
  names(bal) <- c("Control", "Corrupt", "Honest", "Totals")
  bal_table <- t(data.frame(N=bal))
  
  xtable(bal_table) %>%
    print(hline.after = c(),
          only.contents = TRUE,
          file = paste0(tabpath, "treats_table_", study[i], ".tex"))
  
  varnames <- str_replace(covs, ".*_", "")  
  
  desc_df <-   
    combined_trust %>%
    filter(experiment == study[i]) %>%
    select(covs) %>%
    summarise_all(., funs(MEAN = mean(., na.rm = TRUE),
                          SD = sd(., na.rm = TRUE),
                          MEDIAN = median(., na.rm = TRUE),
                          Min = min(., na.rm = TRUE),
                          Max = max(., na.rm = TRUE)
    )) %>%
    gather(., key = "variable") %>%
    mutate(stat = c(rep("Mean", length(covs)),
                    rep("SD", length(covs)), 
                    rep("Median", length(covs)), 
                    rep("Min", length(covs)),
                    rep("Max", length(covs))),
           variable = factor(str_extract(variable, varnames), levels = varnames,
                             labels = c("Female", "Age", "College degree",
                                        "Employed", "Asian", "Black", 
                                        "Hispanic", "White", "Democrat", 
                                        "Republican"))) %>%
    group_by(stat) %>%
    spread(., key = stat, value = value) %>%
    select(variable, Mean, SD, Median, Min, Max) %>%
    mutate(Min = as.integer(Min), 
           Max = as.integer(Max))
  
  xtable(desc_df) %>%
    print(hline.after = c(),
          only.contents = TRUE,
          include.rownames = FALSE, 
          include.colnames = FALSE,
          file = paste0(tabpath, "covs_table_", study[i], ".tex"))
  
}

# Export reliability coefficient for the combined trust measures.
omega_likert <- omega_fedspend <- list()
n <- numeric()
for(i in 1:length(study)){
  omega_likert[[i]] <- 
    combined_trust %>%
    filter(experiment == study[i] & Z == "Placebo") %>% 
    dplyr::select(starts_with("fmptrust_"), -fmptrust_scale, 
                  -fmptrust_index) %>%
    na.omit() %>%
    ci.reliability()
  
  omega_fedspend[[i]] <- 
    combined_trust %>%
    filter(experiment == study[i] & Z == "Placebo") %>% 
    dplyr::select(fedspend_welfare, fedspend_blacks, fedspend_foodstamp,
                  fedspend_homeless) %>%
    na.omit() %>%
    ci.reliability()
  
  n[i] <- 
    combined_trust %>% 
    filter(experiment == study[i] & Z == "Placebo") %>%
    nrow()
}

likert_df <- 
  data.frame(round(bind_rows(omega_likert)[, 1:4], 2), study) %>%
  mutate(n = as.integer(n), 
         study = c("Experiment 1", "Experiment 2", 
                   "Experiment 4", "Experiment 5", 
                   "Experiment 3")) %>%
  select(est, se, ci.lower, ci.upper, n, study) %>% 
  arrange(study)

# Export reliability coefficient for the combined policy measures.
omega_list <- list()
for(i in 1:length(study)){
  omega_list[[i]] <- 
    combined_trust %>%
    filter(experiment == study[i]) %>% 
    select(starts_with("fedspend_"), -fedspend_scale, -fedspend_index) %>%
    na.omit() %>%
    ci.reliability()
}

data.frame(round(bind_rows(omega_list)[, 1:4], 2), study)

fedspend_df <- 
  data.frame(round(bind_rows(omega_fedspend)[, 1:4], 2), study) %>%
  mutate(n = as.integer(n), 
         study = c("Experiment 1", "Experiment 2", 
                   "Experiment 4", "Experiment 5", 
                   "Experiment 3")) %>%
  select(est, se, ci.lower, ci.upper, n, study) %>% 
  arrange(study)

xtable(likert_df) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "omega_table_likert", ".tex"))

xtable(fedspend_df) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "omega_table_fedspend", ".tex"))

#####--------------------------Table S14-----------------------------------#####

## NB: run set.seed() and everything below concurrently for reproducability 
set.seed(123) 


## RI for Experiment 1
exp1_df <- 
  combined_trust %>% 
  filter(!is.na(fedspend_index), experiment == "MTurk_Pols14")

design_exp1 <- declare_ra(N = nrow(exp1_df), num_arms = 3, simple = TRUE,
                          conditions = c("Placebo", "Corrupt", "Honest"))

ri_exp1 <- 
  conduct_ri(
    model_1 = fedspend_index ~ 1, # restricted model
    model_2 = fedspend_index ~ Z, # unrestricted model
    declaration = design_exp1,
    sharp_hypothesis = 0,
    data = exp1_df,
    sims = 5000
  )

## RI for Experiment 2
exp2_df <- 
  combined_trust %>% 
  filter(!is.na(fedspend_index), experiment == "GenPop_Pols14")

block_mat <- as.matrix(table(exp2_df$dem_pid3, exp2_df$Z))
design_exp2 <- declare_ra(N = nrow(exp2_df), num_arms = 3, 
                          blocks = exp2_df$dem_pid3,
                          block_m_each = block_mat,
                          conditions = c("Placebo", "Corrupt", "Honest"))

ri_exp2 <-
  conduct_ri(
    model_1 = fedspend_index ~ 1, # restricted model
    model_2 = fedspend_index ~ Z, # unrestricted model
    declaration = design_exp2,
    sharp_hypothesis = 0,
    data = exp2_df,
    sims = 5000
  )

## RI for Experiment 3
exp3_df <- 
  combined_trust %>% 
  filter(!is.na(fedspend_index), experiment == "MTurk_Pols17")

design_exp3 <- declare_ra(N = nrow(exp3_df), num_arms = 3, simple = TRUE,
                          conditions = c("Placebo", "Corrupt", "Honest"))


ri_exp3 <-
  conduct_ri(
    model_1 = fedspend_index ~ 1, # restricted model
    model_2 = fedspend_index ~ Z, # unrestricted model
    declaration = design_exp3,
    sharp_hypothesis = 0,
    data = exp3_df,
    sims = 5000
  )


## Combine tables:
ri_pooled_df <- 
  bind_rows(ri_exp1$sims_df %>% mutate(experiment = "Experiment 1"), 
            ri_exp2$sims_df %>% mutate(experiment = "Experiment 2"), 
            ri_exp3$sims_df %>% mutate(experiment = "Experiment 3"))

ri_table <- 
  ri_pooled_df %>%
  group_by(experiment) %>%
  summarise(obs = mean(est_obs), 
            ci_lower = quantile(est_sim, probs = 0.025),
            ci_upper = quantile(est_sim, probs = 0.975),
            ri_pvalue = mean(est_sim >= est_obs)
  ) 


xtable(ri_table) %>%
  print(hline.after = c(),
        only.contents = TRUE,
        include.rownames = FALSE, 
        include.colnames = FALSE,
        file = paste0(tabpath, "sharp_ri_redist_table.tex"))


#####--------------------------Tables S15-S18------------------------------#####
summary_df <- 
  combined_trust %>%
  group_by(Z, experiment) %>%
  do(tidy(lm_robust(fmptrust_scale ~ 1, weights = wts, data = .))) %>% 
  select(-term) %>% 
  mutate(measure = "Likert scale") %>%
  rbind(combined_trust %>%
          group_by(Z, experiment) %>% 
          do(tidy(lm_robust(trustgov_scale ~ 1, weights = wts, data = .))) %>% 
          select(-term) %>%  
          mutate(measure = "ANES item")) %>%
  ungroup() %>%
  mutate(Z = factor(Z, levels = c("Corrupt", "Placebo", "Honest"))) %>%
  data.frame

# Tables
group_means <-
  summary_df %>%
  mutate(entry = table_entry(est = estimate, se = std.error, digits = 2)) %>%
  rename(group = experiment) %>%
  select(measure, Z, entry, group) %>%
  spread(Z, entry) %>% 
  mutate(group = factor(group, levels = c("MTurk_Pols14", "GenPop_Pols14",
                                          "MTurk_Pols17", "MTurk_NFL14", 
                                          "MTurk_NFL15"),
                        labels = paste("Experiment", 1:5))) %>% 
  select(group, everything()) %>% 
  arrange(group, desc(measure))


xtable(group_means %>% filter(measure == "ANES item")) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "trust_groupmean_anes.tex"))

xtable(group_means %>% filter(measure != "ANES item")) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "trust_groupmean_likert.tex"))

## Repeat for federal spending DVs
summary_df <- 
  combined_trust %>%
  group_by(Z, experiment) %>%
  do(tidy(lm_robust(fedspend_welfare ~ 1, weights = wts, data = .))) %>% 
  select(-term) %>% 
  mutate(measure = "Welfare") %>%
  rbind(combined_trust %>%
          group_by(Z, experiment) %>% 
          do(tidy(lm_robust(fedspend_foodstamp ~ 1, weights = wts, data = .))) %>% 
          select(-term) %>%  
          mutate(measure = "Foodstamps")) %>% 
  rbind(combined_trust %>%
          group_by(Z, experiment) %>% 
          do(tidy(lm_robust(fedspend_blacks ~ 1, weights = wts, data = .))) %>% 
          select(-term) %>%  
          mutate(measure = "Aid to Blacks")) %>% 
  rbind(combined_trust %>%
          group_by(Z, experiment) %>% 
          do(tidy(lm_robust(fedspend_homeless ~ 1, weights = wts, data = .))) %>% 
          select(-term) %>%  
          mutate(measure = "Aid to Homeless")) %>% 
  rbind(combined_trust %>%
          group_by(Z, experiment) %>% 
          do(tidy(lm_robust(fedspend_scale ~ 1, weights = wts, data = .))) %>% 
          select(-term) %>%  
          mutate(measure = "Redistribution Scale")) %>%
  ungroup() %>%
  mutate(Z = factor(Z, levels = c("Corrupt", "Placebo", "Honest"))) %>%
  data.frame

# Tables
group_means <-
  summary_df %>%
  mutate(entry = table_entry(est = estimate, se = std.error, digits = 2)) %>%
  rename(group = experiment) %>%
  select(measure, Z, entry, group) %>%
  spread(Z, entry) %>% 
  mutate(group = factor(group, levels = c("MTurk_Pols14", "GenPop_Pols14",
                                          "MTurk_Pols17", "MTurk_NFL14", 
                                          "MTurk_NFL15"),
                        labels = paste("Experiment", 1:5)),
         measure = factor(measure, 
                          levels = c("Welfare", "Foodstamps", "Aid to Blacks", 
                                     "Aid to Homeless", "Redistribution Scale"))) %>% 
  select(group, everything()) %>% 
  arrange(group, desc(measure))

xtable(group_means %>% filter(measure != "Redistribution Scale")) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "fedspend_groupmean_all.tex"))

xtable(group_means %>% filter(measure == "Redistribution Scale")) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "fedspend_groupmean_index.tex"))

#####--------------------------Tables S19-S20------------------------------#####
itt_d_fit1 <- 
  lm_robust(D ~ Z + dem_pid3 + dem_age + dem_female + dem_college + dem_work +
              dem_income + dem_white + dem_ideo5 + experiment, weights = wts,
            data = est_df) %>% 
  tidy() %>% 
  filter(term %in% c("ZCorrupt", "ZHonest")) %>% 
  mutate(term = ifelse(term == "ZCorrupt", "Corrupt", "Honest"),
         outcome = "Trust in Government (Likert)",
         S = "Political Corruption") %>% 
  mutate_if(is.numeric, round, 3)

itt_y_fit1 <- 
  lm_robust(Y ~ Z + dem_pid3 + dem_age + dem_female + dem_college + dem_work +
              dem_income + dem_white + dem_ideo5 + experiment, weights = wts,
            data = est_df) %>% 
  tidy() %>% 
  filter(term %in% c("ZCorrupt", "ZHonest")) %>% 
  mutate(term = ifelse(term == "ZCorrupt", "Corrupt", "Honest"),
         outcome = "Support for Redistribution",
         S = "Political Corruption") %>% 
  mutate_if(is.numeric, round, 3)

itt_d_fit2 <- 
  lm_robust(D ~ Z + dem_pid3 + dem_age + dem_female + dem_college + dem_work +
              dem_income + dem_white + dem_ideo5 + experiment, weights = wts,
            data = placebo_df) %>% 
  tidy() %>% 
  filter(term %in% c("ZCorrupt", "ZHonest")) %>% 
  mutate(term = ifelse(term == "ZCorrupt", "Corrupt", "Honest"),
         outcome = "Trust in Government (Likert)",
         S = "Non-Political Corruption") %>% 
  mutate_if(is.numeric, round, 3)

itt_y_fit2 <- 
  lm_robust(Y ~ Z + dem_pid3 + dem_age + dem_female + dem_college + dem_work +
              dem_income + dem_white + dem_ideo5 + experiment, weights = wts,
            data = placebo_df) %>% 
  tidy() %>% 
  filter(term %in% c("ZCorrupt", "ZHonest")) %>% 
  mutate(term = ifelse(term == "ZCorrupt", "Corrupt", "Honest"),
         outcome = "Support for Redistribution",
         S = "Non-Political Corruption") %>% 
  mutate_if(is.numeric, round, 3)

# Estimate ITTd instead using ANES trust measure
itt_d_fit3 <- 
  lm_robust(D1 ~ Z + dem_pid3 + dem_age + dem_female + dem_college + dem_work +
              dem_income + dem_white + dem_ideo5 + experiment, weights = wts,
            data = est_df) %>% 
  tidy() %>% 
  filter(term %in% c("ZCorrupt", "ZHonest")) %>% 
  mutate(term = ifelse(term == "ZCorrupt", "Corrupt", "Honest"),
         outcome = "Trust in Government (ANES item)",
         S = "Political Corruption") %>% 
  mutate_if(is.numeric, round, 3)

itt_d_fit4 <- 
  lm_robust(D1 ~ Z + dem_pid3 + dem_age + dem_female + dem_college + dem_work +
              dem_income + dem_white + dem_ideo5 + experiment, weights = wts,
            data = placebo_df) %>% 
  tidy() %>% 
  filter(term %in% c("ZCorrupt", "ZHonest")) %>% 
  mutate(term = ifelse(term == "ZCorrupt", "Corrupt", "Honest"),
         outcome = "Trust in Government (ANES item)",
         S = "Non-Political Corruption") %>% 
  mutate_if(is.numeric, round, 3)

itt_est_df <- bind_rows(itt_d_fit1, itt_d_fit2, itt_d_fit3, itt_d_fit4,
                        itt_y_fit1, itt_y_fit2)

# Export table for SI: 
pooled_regtable_separate <-
  itt_est_df %>% 
  mutate(entry = make_entry(est = estimate, se = std.error, p = p.value),
         outcome = factor(outcome, levels = c("Trust in Government (Likert)",
                                              "Trust in Government (ANES item)",
                                              "Support for Redistribution"))) %>% 
  select(term, outcome, entry, S) %>%
  spread(term, entry) %>%
  select(S, everything()) %>% 
  arrange(desc(S))

xtable(pooled_regtable_separate) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "pooled_regtable_separate.tex"))


##Repeat above analyses without covariate adjustment
itt_d_fit1 <- 
  lm_robust(D ~ Z + experiment,  weights = wts, data = est_df) %>% 
  tidy() %>% 
  filter(term %in% c("ZCorrupt", "ZHonest")) %>% 
  mutate(term = ifelse(term == "ZCorrupt", "Corrupt", "Honest"),
         outcome = "Trust in Government (Likert)",
         S = "Political Corruption") %>% 
  mutate_if(is.numeric, round, 3)

itt_y_fit1 <- 
  lm_robust(Y ~ Z + experiment, weights = wts, data = est_df) %>% 
  tidy() %>% 
  filter(term %in% c("ZCorrupt", "ZHonest")) %>% 
  mutate(term = ifelse(term == "ZCorrupt", "Corrupt", "Honest"),
         outcome = "Support for Redistribution",
         S = "Political Corruption") %>% 
  mutate_if(is.numeric, round, 3)

itt_d_fit2 <- 
  lm_robust(D ~ Z + experiment, weights = wts, data = placebo_df) %>% 
  tidy() %>% 
  filter(term %in% c("ZCorrupt", "ZHonest")) %>% 
  mutate(term = ifelse(term == "ZCorrupt", "Corrupt", "Honest"),
         outcome = "Trust in Government (Likert)",
         S = "Non-Political Corruption") %>% 
  mutate_if(is.numeric, round, 3)

itt_y_fit2 <- 
  lm_robust(Y ~ Z + experiment, weights = wts, data = placebo_df) %>% 
  tidy() %>% 
  filter(term %in% c("ZCorrupt", "ZHonest")) %>% 
  mutate(term = ifelse(term == "ZCorrupt", "Corrupt", "Honest"),
         outcome = "Support for Redistribution",
         S = "Non-Political Corruption") %>% 
  mutate_if(is.numeric, round, 3)

# Estimate ITTd instead using ANES trust measure
itt_d_fit3 <- 
  lm_robust(D1 ~ Z + experiment, weights = wts, data = est_df) %>% 
  tidy() %>% 
  filter(term %in% c("ZCorrupt", "ZHonest")) %>% 
  mutate(term = ifelse(term == "ZCorrupt", "Corrupt", "Honest"),
         outcome = "Trust in Government (ANES item)",
         S = "Political Corruption") %>% 
  mutate_if(is.numeric, round, 3)

itt_d_fit4 <- 
  lm_robust(D1 ~ Z + experiment, weights = wts, data = placebo_df) %>% 
  tidy() %>% 
  filter(term %in% c("ZCorrupt", "ZHonest")) %>% 
  mutate(term = ifelse(term == "ZCorrupt", "Corrupt", "Honest"),
         outcome = "Trust in Government (ANES item)",
         S = "Non-Political Corruption") %>% 
  mutate_if(is.numeric, round, 3)

itt_est_df <- bind_rows(itt_d_fit1, itt_d_fit2, itt_d_fit3, itt_d_fit4,
                        itt_y_fit1, itt_y_fit2)

# Export table for SI: 
pooled_regtable_separate_unadjusted <-
  itt_est_df %>% 
  mutate(entry = make_entry(est = estimate, se = std.error, p = p.value),
         outcome = factor(outcome, levels = c("Trust in Government (Likert)",
                                              "Trust in Government (ANES item)",
                                              "Support for Redistribution"))) %>% 
  select(term, outcome, entry, S) %>%
  spread(term, entry) %>%
  select(S, everything()) %>% 
  arrange(desc(S)) 

xtable(pooled_regtable_separate_unadjusted) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "pooled_regtable_separate_unadjusted.tex"))


# Compare 2SLS estimates. Create two separate instruments
est_df <- 
  est_df %>% 
  mutate(W1 = ifelse(Z == "Honest", 1, NA),
         W1 = ifelse(Z == "Placebo", 0, W1),
         W2 = ifelse(Z == "Corrupt", 1, NA),
         W2 = ifelse(Z == "Placebo", 0, W2))

late_fit_linear <- 
  iv_robust(Y ~ D + dem_pid3 + dem_age + dem_female + dem_college + dem_income + 
              dem_work + dem_white + dem_ideo5 + experiment | Z1 + dem_pid3 + 
              dem_age + dem_female + dem_college + dem_income + dem_work + 
              dem_white + dem_ideo5 + experiment, weights = wts,
            data = est_df) %>% 
  tidy()

late_fit_honest <- 
  iv_robust(Y ~ D + dem_pid3 + dem_age + dem_female + dem_college + dem_income + 
              dem_work + dem_white + dem_ideo5 + experiment | W1 + dem_pid3 + 
              dem_age + dem_female + dem_college + dem_income + dem_work + 
              dem_white + dem_ideo5 + experiment, weights = wts,
            data = est_df) %>% 
  tidy()

late_fit_corrupt <- 
  iv_robust(Y ~ D + dem_pid3 + dem_age + dem_female + dem_college + dem_income + 
              dem_work + dem_white + dem_ideo5 + experiment | W2 + dem_pid3 + 
              dem_age + dem_female + dem_college + dem_income + dem_work + 
              dem_white + dem_ideo5 + experiment, weights = wts,
            data = est_df) %>% 
  tidy()

#####--------------------------Tables S23-S24------------------------------#####

# Welfare 
itt_y_welfare <- 
  lm_robust(fedspend_welfare ~ Z1 + dem_pid3 + dem_age + dem_female + 
              dem_college + dem_work + dem_income + dem_white + dem_ideo5 +
              experiment, weights = wts, data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "Z1")

late_welfare <- 
  iv_robust(fedspend_welfare ~ D + dem_pid3 + dem_age + dem_female + 
              dem_college + dem_work + dem_income + dem_white + dem_ideo5 +
              experiment | Z1 + dem_pid3 + dem_age + dem_female + dem_college +
              dem_work + dem_income + dem_white + dem_ideo5 + experiment, 
            weights = wts, data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "D")

# Food Stamps
itt_y_foodstamp <- 
  lm_robust(fedspend_foodstamp ~ Z1 + dem_pid3 + dem_age + dem_female + 
              dem_college + dem_work + dem_income + dem_white + dem_ideo5 +
              experiment, weights = wts, data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "Z1")


late_foodstamp <- 
  iv_robust(fedspend_foodstamp ~ D + dem_pid3 + dem_age + dem_female + 
              dem_college + dem_work + dem_income + dem_white + dem_ideo5 +
              experiment | Z1 + dem_pid3 + dem_age + dem_female + dem_college +
              dem_work + dem_income + dem_white + dem_ideo5 + experiment, 
            weights = wts, data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "D")

# Homeless
itt_y_homeless <- 
  lm_robust(fedspend_homeless ~ Z1 + dem_pid3 + dem_age + dem_female + 
              dem_college + dem_work + dem_income + dem_white + dem_ideo5 +
              experiment, weights = wts, data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "Z1")

late_homeless <- 
  iv_robust(fedspend_homeless ~ D + dem_pid3 + dem_age + dem_female + 
              dem_college + dem_work + dem_income + dem_white + dem_ideo5 +
              experiment | Z1 + dem_pid3 + dem_age + dem_female + dem_college +
              dem_work + dem_income + dem_white + dem_ideo5 + experiment, 
            weights = wts, data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "D")

# Aid to blacks
itt_y_blacks <- 
  lm_robust(fedspend_blacks ~ Z1 + dem_pid3 + dem_age + dem_female + 
              dem_college + dem_work + dem_income + dem_white + dem_ideo5 +
              experiment, weights = wts, data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "Z1")

late_blacks <- 
  iv_robust(fedspend_blacks ~ D + dem_pid3 + dem_age + dem_female + 
              dem_college + dem_work + dem_income + dem_white + dem_ideo5 +
              experiment | Z1 + dem_pid3 + dem_age + dem_female + dem_college +
              dem_work + dem_income + dem_white + dem_ideo5 + experiment, 
            weights = wts, data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "D")

# Export table for SI. 
itt_late_sep_table <- 
  bind_rows(itt_y_welfare, itt_y_foodstamp, itt_y_homeless, itt_y_blacks,
            late_welfare, late_foodstamp, late_homeless, late_blacks) %>% 
  mutate(outcome = factor(outcome, levels = c("fedspend_welfare",
                                              "fedspend_foodstamp",
                                              "fedspend_homeless",
                                              "fedspend_blacks"),
                          labels = c("Welfare", "Food Stamps", 
                                     "Aid to Homeless", "Aid to Blacks")),
         term = factor(ifelse(term == "Z1", "Reduced Form", "2SLS"), 
                       levels = c("Reduced Form", "2SLS")),
         entry = make_entry(est = estimate, se = std.error, p = p.value)) %>% 
  select(term, entry, outcome) %>% 
  spread(outcome, entry)


xtable(itt_late_sep_table)  %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "itt_late_sep_table.tex"))


# Repeat without covariate adjustment

# Welfare 
itt_y_welfare <- 
  lm_robust(fedspend_welfare ~ Z1  + experiment, weights = wts,
            data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "Z1")

late_welfare <- 
  iv_robust(fedspend_welfare ~ D + experiment | Z1 + experiment, weights = wts, 
            data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "D")

# Food Stamps
itt_y_foodstamp <- 
  lm_robust(fedspend_foodstamp ~ Z1 + experiment, weights = wts,
            data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "Z1")


late_foodstamp <- 
  iv_robust(fedspend_foodstamp ~ D  + experiment | Z1 + experiment, 
            weights = wts, data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "D")

# Homeless
itt_y_homeless <- 
  lm_robust(fedspend_homeless ~ Z1 + experiment, weights = wts,
            data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "Z1")

late_homeless <- 
  iv_robust(fedspend_homeless ~ D + experiment | Z1 + experiment, weights = wts,
            data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "D")

# Aid to blacks
itt_y_blacks <- 
  lm_robust(fedspend_blacks ~ Z1 + experiment, weights = wts,
            data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "Z1")

late_blacks <- 
  iv_robust(fedspend_blacks ~ D  + experiment | Z1 + experiment, weights = wts,
            data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term == "D")

# Export table for SI. 
itt_late_sep_table <- 
  bind_rows(itt_y_welfare, itt_y_foodstamp, itt_y_homeless, itt_y_blacks,
            late_welfare, late_foodstamp, late_homeless, late_blacks) %>% 
  mutate(outcome = factor(outcome, levels = c("fedspend_welfare",
                                              "fedspend_foodstamp",
                                              "fedspend_homeless",
                                              "fedspend_blacks"),
                          labels = c("Welfare", "Food Stamps", 
                                     "Aid to Homeless", "Aid to Blacks")),
         term = factor(ifelse(term == "Z1", "Reduced Form", "2SLS"), 
                       levels = c("Reduced Form", "2SLS")),
         entry = make_entry(est = estimate, se = std.error, p = p.value)) %>% 
  select(term, entry, outcome) %>% 
  spread(outcome, entry)


xtable(itt_late_sep_table)  %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "itt_late_sep_table_unadjusted.tex"))

#####--------------------------Table S24-----------------------------------#####

# ITTd:
itt_d_fit1 <- 
  lm_robust(D ~ Z1 + experiment, weights = wts, data = est_df) 

itt_d_fit2 <- 
  lm_robust(D ~ Z1 + experiment, weights = wts, data = placebo_df) 

# ITTy:
itt_y_fit1 <- 
  lm_robust(Y ~ Z1 + experiment, weights = wts, data = est_df) 

itt_y_fit2 <- 
  lm_robust(Y ~ Z1 + experiment, weights = wts, data = placebo_df) 

itt_d_est <- 
  itt_d_fit1 %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c('(Intercept)', "Z1")) %>% 
  mutate(S = "Political Corruption") %>% 
  bind_rows(
    itt_d_fit2 %>% 
      tidy() %>% 
      mutate_if(is.numeric, round, 3) %>% 
      filter(term %in% c('(Intercept)', "Z1")) %>% 
      mutate(S = "Non-Political Corruption")
  ) %>% 
  mutate(estimand = "Trust in Government (First Stage)")


itt_y_est <-
  itt_y_fit1 %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c('(Intercept)', "Z1")) %>% 
  mutate(S = "Political Corruption") %>% 
  bind_rows(
    itt_y_fit2 %>% 
      tidy() %>% 
      mutate_if(is.numeric, round, 3) %>% 
      filter(term %in% c('(Intercept)', "Z1")) %>% 
      mutate(S = "Non-Political Corruption")
  ) %>% 
  mutate(estimand = "Support for Redistribution (Reduced Form)")


itt_est_df <- bind_rows(itt_d_est, itt_y_est)

# Estimate LATE of D on Y in politics experiments
late_fit1 <- 
  iv_robust(Y ~ D + experiment | Z1 + experiment, weights = wts,
            data = est_df) 

late_est_df <- 
  late_fit1 %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c('(Intercept)', "D")) %>% 
  mutate(estimand = "LATE", S = "Political Corruption") 

# Export table for SI: 
pooled_regtable <-
  bind_rows(itt_est_df, late_est_df) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(entry = make_entry(est = estimate, se = std.error, p = p.value)) %>% 
  select(estimand, entry, S) %>%
  spread(estimand, entry) %>%
  arrange(desc(S)) %>% 
  select(S, `Trust in Government (First Stage)`, 
         `Support for Redistribution (Reduced Form)`, LATE)

pooled_regtable[2,4] <- "-"

xtable(pooled_regtable) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "pooled_regtable_unadjusted.tex"))

#####--------------------------Table S25-----------------------------------#####

# Test of instrument strength: 
summary(lm(D ~ Z1, weights = wts, data = est_df)) # Likert trust
summary(lm(D1 ~ Z1, weights = wts, data = est_df)) # ANES trust

# Compare to placebo experiments
summary(lm(D ~ Z1, data = placebo_df)) # Likert trust
summary(lm(D1 ~ Z1, data = placebo_df)) # Likert trust # ANES trust

## Loop through DVs, estimating ITT, LATE 
covariates <- c("dem_pid3", "dem_age", "dem_female", "dem_college", 
                "dem_income", "dem_white", "dem_ideo5", "dem_work")

all_dvs <- c("fmptrust_scale", "trustgov_scale", "fedspend_scale")

S <- c("MTurk_Pols14", "GenPop_Pols14", "MTurk_Pols17", "MTurk_NFL14", 
       "MTurk_NFL15")

gg_itt <- list()
gg_itt_df <- list()

for(j in 1:length(S)){
  for(i in 1:length(all_dvs)) {
    
    Y <- all_dvs[i]
    X <- covariates
    
    gg_itt[[i]] <- 
      lm_robust(as.formula(paste(Y, "~ Z1 +", paste(X, collapse = "+"))), 
                weights = wts, data = combined_trust %>% 
                  filter(experiment == S[j])) %>%
      tidy() %>%
      mutate_if(is.numeric, round, 3) %>%
      filter(term %in% c("(Intercept)", "Z1")) %>%
      mutate(estimator = "OLS",
             term = ifelse(term == "(Intercept)", "Control", "Treatment"))
    
  }
  gg_itt_df[[j]] <- bind_rows(gg_itt) %>%  mutate(experiment = paste(S[j]))
}


all_itt_df <- 
  bind_rows(gg_itt_df) %>% 
  filter(term == "Treatment") %>% 
  mutate(study = NA, 
         study = ifelse(experiment == "MTurk_Pols14", "Experiment 1", 
                        study),
         study = ifelse(experiment == "GenPop_Pols14", "Experiment 2", 
                        study),
         study = ifelse(experiment == "MTurk_NFL14", "Experiment 4", 
                        study),
         study = ifelse(experiment == "MTurk_NFL15", "Experiment 5", 
                        study),
         study = ifelse(experiment == "MTurk_Pols17", "Experiment 3", 
                        study),
         study = factor(study, levels = c("Experiment 1", 
                                          "Experiment 2", 
                                          "Experiment 3",
                                          "Experiment 4", 
                                          "Experiment 5")),
         estimator = factor(outcome, levels = c("trustgov_scale",
                                                "fmptrust_scale", 
                                                "fedspend_scale"),
                            labels = c("First Stage (ANES item)",
                                       "First Stage (Likert Scale)",
                                       "Reduced Form (Redistribution Index)")))


S <- c("MTurk_Pols14", "GenPop_Pols14", "MTurk_Pols17")
gg_late_fmp <- list()
gg_late_anes <- list()

for(j in 1:length(S)){
  X <- covariates
  
  gg_late_fmp[[j]] <- 
    iv_robust(as.formula(paste("fedspend_scale ~ fmptrust_scale + ", 
                               paste(X, collapse = "+"), "| Z1 + ", 
                               paste(X, collapse = "+"))),
              weights = wts,
              data = combined_trust %>% filter(experiment == S[j])) %>%
    tidy() %>%
    mutate_if(is.numeric, round, 3) %>%
    filter(term %in% c("(Intercept)", "fmptrust_scale")) %>%
    mutate(estimator = "2SLS (Likert Scale)", 
           term = ifelse(term == "(Intercept)", "Control", "Treatment"),
           experiment = paste(S[j]))
  
  gg_late_anes[[j]] <- 
    iv_robust(as.formula(paste("fedspend_scale ~ trustgov_scale + ", 
                               paste(X, collapse = "+"),
                               "| Z1 + ", paste(X, collapse = "+"))),
              weights = wts,
              data = combined_trust %>% filter(experiment == S[j])) %>%
    tidy() %>%
    mutate_if(is.numeric, round, 3) %>%
    filter(term %in% c("(Intercept)", "trustgov_scale")) %>%
    mutate(estimator = "2SLS (ANES item)", 
           term = ifelse(term == "(Intercept)", "Control", "Treatment"),
           experiment = paste(S[j]))
  
}

all_late_df <- 
  bind_rows(gg_late_fmp, gg_late_anes) %>% 
  filter(term == "Treatment") %>% 
  mutate(study = NA, 
         study = ifelse(experiment == "MTurk_Pols14", "Experiment 1", 
                        study),
         study = ifelse(experiment == "GenPop_Pols14", "Experiment 2", 
                        study),
         study = ifelse(experiment == "MTurk_NFL14", "Experiment 4", 
                        study),
         study = ifelse(experiment == "MTurk_NFL15", "Experiment 5", 
                        study),
         study = ifelse(experiment == "MTurk_Pols17", "Experiment 3", 
                        study),
         study = factor(study, levels = c("Experiment 1", 
                                          "Experiment 2", 
                                          "Experiment 3",
                                          "Experiment 4", 
                                          "Experiment 5")))

all_itt_late_table <- 
  bind_rows(all_itt_df, all_late_df) %>% 
  mutate(entry = make_entry(est = estimate, se = std.error, p = p.value),
         estimator = factor(estimator, levels = c("First Stage (ANES item)",
                                                  "First Stage (Likert Scale)",
                                                  "Reduced Form (Redistribution Index)",
                                                  "2SLS (ANES item)", 
                                                  "2SLS (Likert Scale)"))) %>% 
  dplyr::select(entry, estimator, study) %>%
  spread(estimator, entry) 


all_itt_late_table[4:5, 5:6] <- "-"
xtable(all_itt_late_table) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "all_itt_late_table.tex"))


#####--------------------------Tables S26-S29------------------------------#####

# Descriptives reported in paper
summary(combined_trust$dem_ideo5) # mean/median ideology score

# Diffs on support for redistribuiton for employed v unemployed folks
combined_trust %>% 
  filter(Z == "Placebo") %>% 
  mutate(dem_work = as.logical(dem_work)) %>% 
  group_by(dem_work) %>% 
  summarise(avg = mean(fedspend_scale, na.rm = T))


# Race: White v. Non-White
het_late_race <- 
  iv_robust(Y ~ D*dem_white + experiment | Z1*dem_white + experiment,
            weights = wts, data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c("(Intercept)", "D", "dem_white", "D:dem_white")) %>% 
  mutate(outcome = "Support for Redistribution",
         model = "Race",
         estimator = "2SLS")

# First stage
het_itt_d_race <- lm_robust(D ~ Z1*dem_white + experiment, weights = wts,
                            data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c("(Intercept)", "Z1", "dem_white", "Z1:dem_white")) %>% 
  mutate(outcome = "Trust in Government",
         model = "Race",
         estimator = "First Stage")

# Reduced form
het_itt_y_race <- lm_robust(Y ~ Z1*dem_white + experiment, weights = wts,
                            data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c("(Intercept)", "Z1", "dem_white", "Z1:dem_white")) %>% 
  mutate(outcome = "Support for Redistribution",
         model = "Race",
         estimator = "Reduced Form")

# Employed v. Unenmployed
het_late_emp <- 
  iv_robust(Y ~ D*dem_work + experiment | Z1*dem_work + experiment, 
            weights = wts, data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c("(Intercept)", "D", "dem_work", "D:dem_work")) %>% 
  mutate(outcome = "Support for Redistribution",
         model = "Work",
         estimator = "2SLS")

# First stage
het_itt_d_emp <- lm_robust(D ~ Z1*dem_work + experiment, weights = wts,
                           data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c("(Intercept)", "Z1", "dem_work", "Z1:dem_work")) %>% 
  mutate(outcome = "Trust in Government",
         model = "Work",
         estimator = "First Stage")

# Reduced form
het_itt_y_emp <- lm_robust(Y ~ Z1*dem_work + experiment, weights = wts,
                           data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c("(Intercept)", "Z1", "dem_work", "Z1:dem_work")) %>% 
  mutate(outcome = "Support for Redistribution",
         model = "Work",
         estimator = "Reduced Form")


# Ideology
het_late_ideo <- 
  iv_robust(Y ~ D*dem_ideo5 + experiment | Z1*dem_ideo5 + experiment, 
            weights = wts, data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c("(Intercept)", "D", "dem_ideo5", "D:dem_ideo5")) %>% 
  mutate(outcome = "Support for Redistribution",
         model = "Ideology",
         estimator = "2SLS")

# First stage
het_itt_d_ideo <- lm_robust(D ~ Z1*dem_ideo5 + experiment, weights = wts,
                            data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c("(Intercept)", "Z1", "dem_ideo5", "Z1:dem_ideo5")) %>% 
  mutate(outcome = "Trust in Government",
         model = "Ideology",
         estimator = "First Stage")

# Reduced form
het_itt_y_ideo <- lm_robust(Y ~ Z1*dem_ideo5 + experiment, weights = wts,
                            data = est_df) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c("(Intercept)", "Z1", "dem_ideo5", "Z1:dem_ideo5")) %>% 
  mutate(outcome = "Support for Redistribution",
         model = "Ideology",
         estimator = "Reduced Form")


# Partisanship
het_late_pid <- 
  iv_robust(Y ~ D*dem_pid3 + experiment | Z1*dem_pid3 + experiment, 
            weights = wts, data = est_df %>% 
              mutate(dem_pid3 = factor(dem_pid3, levels = c("Democrat", 
                                                            "Independent",
                                                            "Republican")))) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c("(Intercept)",  "dem_pid3", "D", "D:dem_pid3Independent",
                     "D:dem_pid3Republican")) %>% 
  mutate(outcome = "Support for Redistribution",
         model = "Partisanship",
         estimator = "2SLS")


# First stage
het_itt_d_pid <- 
  lm_robust(D ~ Z1*dem_pid3 + experiment, weights = wts, data = est_df %>% 
              mutate(dem_pid3 = factor(dem_pid3, levels = c("Democrat",
                                                            "Independent",
                                                            "Republican")))) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c("(Intercept)",  "dem_pid3", "Z1", "Z1:dem_pid3Independent",
                     "Z1:dem_pid3Republican")) %>% 
  mutate(outcome = "Trust in Government",
         model = "Partisanship",
         estimator = "First Stage")

# Reduced form
het_itt_y_pid <- 
  lm_robust(Y ~ Z1*dem_pid3 + experiment, weights = wts,
            data = est_df %>% 
              mutate(dem_pid3 = factor(dem_pid3, levels = c("Democrat",
                                                            "Independent",
                                                            "Republican")))) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c("(Intercept)",  "dem_pid3", "Z1", "Z1:dem_pid3Independent",
                     "Z1:dem_pid3Republican")) %>% 
  mutate(outcome = "Support for Redistribution",
         model = "Partisanship",
         estimator = "Reduced Form")

het_itt_d_df <- bind_rows(het_itt_d_race, het_itt_d_emp, het_itt_d_ideo,
                          het_itt_d_pid)

het_itt_y_df <- bind_rows(het_itt_y_race, het_itt_y_emp, het_itt_y_ideo,
                          het_itt_y_pid)

het_late_df <- bind_rows(het_late_race, het_late_emp, het_late_ideo,
                         het_late_pid)


# Export tables of results. 
het_itt_d_table <- 
  het_itt_d_df %>% 
  filter(term %in% c("Z1", "Z1:dem_white", "Z1:dem_work", "Z1:dem_ideo5",
                     "Z1:dem_pid3Independent", "Z1:dem_pid3Republican")) %>% 
  mutate(entry = make_entry(est = estimate, se = std.error, p = p.value),
         model = factor(model, levels = c("Race", "Work", "Partisanship",
                                          "Ideology")),
         term = factor(term, levels = c("Z1", "Z1:dem_white", "Z1:dem_work", 
                                        "Z1:dem_pid3Independent", 
                                        "Z1:dem_pid3Republican",
                                        "Z1:dem_ideo5"),
                       labels = c("Treatment", "Treatment x White", 
                                  "Treatment x Employed", 
                                  "Treatment x Independent", 
                                  "Treatment x Republican", 
                                  "Treatment x Ideology"))) %>% 
  select(term, entry, model) %>%
  spread(model, entry) 

het_itt_d_table[is.na(het_itt_d_table)] <- "-"


xtable(het_itt_d_table) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "het_itt_d_table.tex"))


het_itt_y_table <- 
  het_itt_y_df %>% 
  filter(term %in% c("Z1", "Z1:dem_white", "Z1:dem_work", "Z1:dem_ideo5",
                     "Z1:dem_pid3Independent", "Z1:dem_pid3Republican")) %>% 
  mutate(entry = make_entry(est = estimate, se = std.error, p = p.value),
         model = factor(model, levels = c("Race", "Work", "Partisanship",
                                          "Ideology")),
         term = factor(term, levels = c("Z1", "Z1:dem_white", "Z1:dem_work", 
                                        "Z1:dem_pid3Independent", 
                                        "Z1:dem_pid3Republican",
                                        "Z1:dem_ideo5"),
                       labels = c("Treatment", "Treatment x White", 
                                  "Treatment x Employed", 
                                  "Treatment x Independent", 
                                  "Treatment x Republican", 
                                  "Treatment x Ideology"))) %>% 
  select(term, entry, model) %>%
  spread(model, entry) 

het_itt_y_table[is.na(het_itt_y_table)] <- "-"

xtable(het_itt_y_table) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "het_itt_y_table.tex"))


het_late_table <- 
  het_late_df %>% 
  filter(term %in% c("D", "D:dem_white", "D:dem_work", "D:dem_ideo5",
                     "D:dem_pid3Independent", "D:dem_pid3Republican")) %>% 
  mutate(entry = make_entry(est = estimate, se = std.error, p = p.value),
         model = factor(model, levels = c("Race", "Work", "Partisanship",
                                          "Ideology")),
         term = factor(term, levels = c("D", "D:dem_white", "D:dem_work", 
                                        "D:dem_pid3Independent", 
                                        "D:dem_pid3Republican",
                                        "D:dem_ideo5"),
                       labels = c("Trust in Government", "Trust x White", 
                                  "Trust x Employed", "Trust x Independent",
                                  "Trust x Republican", "Trust x Ideology"))) %>% 
  select(term, entry, model) %>%
  spread(model, entry) 

het_late_table[is.na(het_late_table)] <- "-"

xtable(het_late_table) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "het_late_table.tex"))

# Reduced form CATEs for unemployed v. employed
lm_robust(Y ~ Z1 + experiment, weights = wts, 
          data = est_df %>% filter(dem_work == 0))  %>% 
  tidy() %>% 
  filter(term == "Z1") %>% 
  mutate_if(is.numeric, round, 3)

lm_robust(Y ~ Z1 + experiment, weights = wts,
          data = est_df %>% filter(dem_work == 1))  %>% 
  tidy() %>% 
  filter(term == "Z1") %>% 
  mutate_if(is.numeric, round, 3)


iv_robust(Y ~ D + experiment | Z1 + experiment, 
          weights = wts, data = est_df %>% filter(dem_work == 0)) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c("(Intercept)", "D", "dem_work", "D:dem_work")) %>% 
  mutate(outcome = "Support for Redistribution",
         model = "Work",
         estimator = "2SLS")

iv_robust(Y ~ D + experiment | Z1 + experiment, 
          weights = wts, data = est_df %>% filter(dem_work == 1)) %>% 
  tidy() %>% 
  mutate_if(is.numeric, round, 3) %>% 
  filter(term %in% c("(Intercept)", "D", "dem_work", "D:dem_work")) %>% 
  mutate(outcome = "Support for Redistribution",
         model = "Work",
         estimator = "2SLS")

## Adjusting for multiple comparisons
comparison_df <- 
  bind_rows(het_itt_d_df, het_itt_y_df, het_late_df) %>% 
  filter(term %in% c("Z1:dem_white", "Z1:dem_work", "Z1:dem_pid3Independent",
                     "Z1:dem_pid3Republican", "Z1:dem_ideo5", "D:dem_white", 
                     "D:dem_work", "D:dem_pid3Independent",
                     "D:dem_pid3Republican", "D:dem_ideo5")) %>% 
  mutate(p.bh = p.adjust(p.value, method = "BH"),
         model = factor(model, levels = c("Race", "Work", "Partisanship",
                                          "Ideology")),
         term = factor(term, levels = c("Z1:dem_white", "Z1:dem_work", 
                                        "Z1:dem_pid3Independent", 
                                        "Z1:dem_pid3Republican",
                                        "Z1:dem_ideo5", "D:dem_white", 
                                        "D:dem_work", "D:dem_pid3Independent", 
                                        "D:dem_pid3Republican", "D:dem_ideo5"),
                       labels = c("Treatment x White", "Treatment x Employed", 
                                  "Treatment x Independent", 
                                  "Treatment x Republican", 
                                  "Treatment x Ideology", "Trust x White", 
                                  "Trust x Employed", "Trust x Independent",
                                  "Trust x Republican", "Trust x Ideology"))) %>% 
  select(term, model, estimator, p.value, p.bh)

xtable(comparison_df %>% filter(estimator == "First Stage")) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "first_stage_comparison_table.tex"))

xtable(comparison_df %>% filter(estimator == "Reduced Form")) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "reduced_form_comparison_table.tex"))

xtable(comparison_df %>% filter(estimator == "2SLS")) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "2sls_comparison_table.tex"))

# Model comparisons:
est_df$dem_independent <- ifelse(est_df$dem_pid3 == "Independent", 1, 0)
est_df$experiment <- factor(est_df$experiment, levels = c("GenPop_Pols14",
                                                          "MTurk_Pols14",
                                                          "MTurk_Pols17"))

# First Stage
first_stage_model1 <- 
  lm_robust(D ~ Z1 + dem_white + dem_work + dem_ideo5 + dem_independent +
              dem_republican + experiment, weights = wts, data = est_df)

first_stage_model2 <- 
  lm_robust(D ~ Z1*dem_white + Z1*dem_work + Z1*dem_ideo5 + 
              Z1*dem_independent + Z1*dem_republican + experiment, 
            weights = wts, data = est_df)


(first_stage_test <- lmtest::waldtest(first_stage_model1, first_stage_model2, 
                                      test = "F"))

# Reduced form
reduced_form_model1 <- 
  lm_robust(Y ~ Z1 + dem_white + dem_work + dem_ideo5 + dem_independent +
              dem_republican + experiment, weights = wts, data = est_df)

reduced_form_model2 <- 
  lm_robust(Y ~  Z1*dem_white + Z1*dem_work + Z1*dem_ideo5 + 
              Z1*dem_independent + Z1*dem_republican + experiment, 
            weights = wts, data = est_df)


(reduced_form_test <- lmtest::waldtest(reduced_form_model1, reduced_form_model2, 
                                       test = "F"))



#####--------------------------Tables S30-S31------------------------------#####
covariates <- c("dem_pid3", "dem_age", "dem_female", "dem_college", 
                "dem_income", "dem_white", "dem_ideo5", "dem_work")

all_dvs <- c("govrole_services", "govrole_jobs", "govrole_blacks",
             "govrole_healthcare", "govspend_defense", "fedspend_foreignaid",
             "fedspend_highways", "fedspend_ss", "fedspend_enviro", 
             "fedspend_school", "fedspend_crime", "fedspend_immig")

other_df <- 
  combined_trust %>% 
  filter(study_type == "Politics") %>% 
  mutate_at(all_dvs, standardize)

gg_itt <- list()

for(i in 1:length(all_dvs)) {
  
  Y <- all_dvs[i]
  X <- covariates
  
  gg_itt[[i]] <- 
    lm_robust(as.formula(paste(Y, "~ Z1 +", paste(X, collapse = "+"))), 
              weights = wts, data = other_df) %>%
    tidy() %>%
    mutate_if(is.numeric, round, 3) %>%
    filter(term %in% c("(Intercept)", "Z1")) %>%
    mutate(estimator = "OLS",
           term = ifelse(term == "(Intercept)", "Control", "Treatment"))
  
}


other_itt_df <- 
  bind_rows(gg_itt) %>% 
  filter(term == "Treatment")

other_n <- 
  other_df %>% 
  select(all_dvs) %>% 
  summarise_all(function(x) paste(3837 - sum(is.na(x))))

other_mean <- 
  other_df %>% 
  select(all_dvs) %>% 
  summarise_all(function(x) sprintf("%.2f", mean(x, na.rm = T))) 

other_itt_table <- 
  other_itt_df %>% 
  mutate(entry = make_entry(est = estimate, se = std.error, p = p.value)) %>% 
  dplyr::select(entry, outcome) %>%
  spread(outcome, entry) %>% 
  bind_rows(other_mean, other_n)

rownames(other_itt_table) <- c("Treatment", "Corrupt mean", "Observations")


# Table S30
xtable(other_itt_table %>% 
         select(fedspend_foreignaid, starts_with("govrole_"), 
                govspend_defense)) %>%
  print(include.rownames = TRUE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "govrole_table.tex"))

# Table S31
xtable(other_itt_table %>% 
         select(fedspend_ss, fedspend_enviro, fedspend_crime,
                fedspend_highways, fedspend_school, fedspend_immig)) %>%
  print(include.rownames = TRUE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "fedspend_distributive_table.tex"))

#####--------------------------Tables S32-S34------------------------------#####
covariates <- c("dem_pid3", "dem_age", "dem_female", "dem_college", 
                "dem_income", "dem_white", "dem_ideo5", "dem_work")

all_dvs <- c("bintrust_family", "bintrust_friends", "bintrust_scientists",
             "bintrust_neighbors", "bintrust_unis", "bintrust_media",
             "bintrust_stranger", "bintrust_police", "bintrust_fedgov", 
             "bintrust_govadmins","bintrust_politicians", "bintrust_stategov",
             "bintrust_localgov")

gg_itt <- list()

for(i in 1:length(all_dvs)) {
  
  Y <- all_dvs[i]
  X <- covariates
  
  gg_itt[[i]] <- 
    lm_robust(as.formula(paste(Y, "~ Z1 +", paste(X, collapse = "+"))), 
              weights = wts, data = combined_trust %>% 
                filter(study_type == "Politics")) %>%
    tidy() %>%
    mutate_if(is.numeric, round, 3) %>%
    filter(term %in% c("(Intercept)", "Z1")) %>%
    mutate(estimator = "OLS",
           term = ifelse(term == "(Intercept)", "Control", "Treatment"))
  
}


bintrust_itt_df <- 
  bind_rows(gg_itt) %>% 
  filter(term == "Treatment")

bintrust_n <- 
  combined_trust %>% 
  filter(study_type == "Politics", experiment != "MTurk_Pols17") %>% 
  select(all_dvs) %>% 
  summarise_all(function(x) paste(1976-sum(is.na(x))))

bintrust_mean <- 
  combined_trust %>% 
  filter(study_type == "Politics", experiment != "MTurk_Pols17") %>% 
  select(all_dvs) %>% 
  summarise_all(function(x) sprintf("%.2f", mean(x, na.rm = T))) 

bintrust_itt_table <- 
  bintrust_itt_df %>% 
  mutate(entry = make_entry(est = estimate, se = std.error, p = p.value)) %>% 
  dplyr::select(entry, outcome) %>%
  spread(outcome, entry) %>% 
  bind_rows(bintrust_mean, bintrust_n)

rownames(bintrust_itt_table) <- c("Treatment", "Corrupt mean", "Observations")

xtable(bintrust_itt_table %>% select(bintrust_family, bintrust_friends, 
                                     bintrust_neighbors, bintrust_stranger)) %>%
  print(include.rownames = TRUE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "bintrust_social_table.tex"))

xtable(bintrust_itt_table %>% select(bintrust_media, bintrust_police, 
                                     bintrust_scientists, bintrust_unis)) %>%
  print(include.rownames = TRUE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "bintrust_orgs_table.tex"))

xtable(bintrust_itt_table %>% select(bintrust_fedgov, bintrust_govadmins, 
                                     bintrust_politicians, bintrust_stategov,
                                     bintrust_localgov)) %>%
  print(include.rownames = TRUE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "bintrust_gov_table.tex"))

#####--------------------------Tables S35-S38------------------------------#####

## Read in datase from AER paper by Kuziemko et al. 
## (Link: https://www.aeaweb.org/articles?id=10.1257/aer.20130360)
## The recoding below is the R translation of what's in ineq_Trust.do
aer_df <- 
  haven::read_dta("~/Dropbox/pipeline/CombinedTrust/Archive/Kuziemko/2013-0360_data/Inequality_AER_Trust.dta") %>% 
  filter(age >= 18, endofsurvey != "") 
  # Their Stata code gets rid of 6 folks who are under 18 & 545 folks
  # who do not finish survey 

## Create binary treatment indicator:
aer_df$treat1 <- NA
aer_df$treat1[aer_df$group == "Negative Trust Treatment"] <- 1
aer_df$treat1[aer_df$group == "Control"] <- 0

table(aer_df$treat1)


## Make it binary 
aer_df$trust_gov_yes <- NA
aer_df$trust_gov_yes <- as.numeric(aer_df$trust_government <= 2)

## 88 people in sample "trust" the government. 
table(aer_df$trust_gov_yes)

## Examine instrument strength (f-stat from 1st stage regression)
summary(lm(trust_gov_yes ~ treat1, data = aer_df))
              
## NB: the reported point estimate for first stage effect (Table 8, Col 1) is a 
## bit different. It looks like it's from a covariate adjusted regression w/ 
## state fixed-effects. One of the covariates, HH income is missing 2 
## observations, which could also explain why the reported # of observations is 
## 899 rather than 901

## Can exactly match the control group mean reported in Table 8 if these 2 
## units are removed: 
with(aer_df, mean(trust_gov_yes[!is.na(hhincome) & treat1 == 0]))

## These 2 folks did not "trust" the government 
with(aer_df, table(is.na(hhincome), trust_gov_yes))

## Flip around government trust so it's in increasing direction
aer_df$trust_gov <- NA
aer_df$trust_gov <- 5 - aer_df$trust_government

## First-stage effect when DV is not recoded to binary indicator
summary(lm(trust_gov ~ treat1, data = aer_df))


## Now subset to experiment 3 data and code Kuziemko et al variables to match
## the coding scheme from the AER paper. 
kuz_df <- 
  combined_trust %>%
  mutate(kuzineq_taxmilnrs = standardize(kuzineq_taxmilnrs),
         kuzpol_minwage = standardize(kuzpol_minwage),
         kuzpol_atp = standardize(kuzpol_atp),
         kuzpol_foodstamp = standardize(kuzpol_foodstamp),
         kuzpol_housepoor = standardize(kuzpol_housepoor),
         kuzineq_ineqprob = standardize(kuzineq_ineqprob),
         kuzineq_povprob = standardize(kuzineq_povprob),
         kuzineq_richdeserve = standardize(kuzineq_richdeserve),
         kuzrank_charity = standardize(kuzrank_charity),
         kuzrank_edu = standardize(kuzrank_edu),
         kuzrank_regulate = standardize(kuzrank_regulate),
         kuzrank_transfer = standardize(kuzrank_transfer),
         kuzrank_taxes = standardize(kuzrank_taxes),
         kuzineq_rog = standardize(kuzineq_rog),
         kuztax_top1 = standardize(kuztax_top1),
         kuztax_next9 = standardize(kuztax_next9),
         kuztax_next40 = standardize(kuztax_next40),
         kuztax_btm50 = standardize(kuztax_btm50),
         kuzpol_housepoor = standardize(kuzpol_housepoor),
         kuzpol_foodstamp = standardize(kuzpol_foodstamp),
         kuzpol_atp = standardize(kuzpol_atp),
         trustgov = as.numeric(trustgov_trustgstd >= 3)) %>%
  filter(experiment == "MTurk_Pols17")

lm_robust(trustgov ~ Z, kuz_df)


covariates <- c("dem_pid3", "dem_age", "dem_female", "dem_college", 
                "dem_income", "dem_white", "dem_ideo5", "dem_work")

all_dvs <- c("kuzineq_taxmilnrs", "kuzpol_minwage", "kuzineq_ineqprob", 
             "kuzineq_rog","kuzineq_povprob", "kuzineq_richdeserve", 
             "kuzrank_charity", "kuzrank_edu", "kuzrank_regulate", 
             "kuzrank_transfer", "kuzrank_taxes", "kuztax_top1", 
             "kuztax_next9", "kuztax_next40", "kuztax_btm50",
             "kuzpol_housepoor", "kuzpol_foodstamp", "kuzpol_atp")

gg_itt <- list()

for(i in 1:length(all_dvs)) {
  
  Y <- all_dvs[i]
  X <- covariates
  
  gg_itt[[i]] <- 
    lm_robust(as.formula(paste(Y, "~ Z1 +", paste(X, collapse = "+"))), 
              data = kuz_df) %>%
    tidy() %>%
    mutate_if(is.numeric, round, 3) %>%
    filter(term %in% c("(Intercept)", "Z1")) %>%
    mutate(estimator = "OLS",
           term = ifelse(term == "(Intercept)", "Control", "Treatment"))
  
}

kuz_itt_df <- 
  bind_rows(gg_itt) %>% 
  filter(term == "Treatment")

kuz_n <- 
  kuz_df %>% 
  select(all_dvs) %>% 
  summarise_all(function(x) paste(1870-sum(is.na(x))))

kuz_mean <- 
  kuz_df %>% 
  select(all_dvs) %>% 
  summarise_all(function(x) sprintf("%.2f", mean(x, na.rm = T))) 

kuz_itt_table <- 
  kuz_itt_df %>% 
  mutate(entry = make_entry(est = estimate, se = std.error, p = p.value)) %>% 
  dplyr::select(entry, outcome) %>%
  spread(outcome, entry) %>% 
  bind_rows(kuz_mean, kuz_n)


rownames(kuz_itt_table) <- c("Treatment Effect", "Corrupt mean", 
                             "Observations")


xtable(kuz_itt_table  %>% select(kuzpol_minwage,
                                 kuzpol_housepoor,
                                 kuzpol_foodstamp,
                                 kuzpol_atp,
                                 kuzineq_rog)) %>%
  print(include.rownames = TRUE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "kuz_redist_table.tex"))


xtable(kuz_itt_table  %>% select(kuzineq_taxmilnrs,
                                 kuzineq_richdeserve,
                                 kuzineq_ineqprob,
                                 kuzineq_povprob)) %>%
  print(include.rownames = TRUE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "kuz_ineq_table.tex"))


xtable(kuz_itt_table  %>% select(kuztax_top1,
                                 kuztax_next9,
                                 kuztax_next40,
                                 kuztax_btm50)) %>%
  print(include.rownames = TRUE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "kuz_tax_table.tex"))


xtable(kuz_itt_table  %>% select(kuzrank_charity,
                                 kuzrank_edu,
                                 kuzrank_transfer,
                                 kuzrank_regulate,
                                 kuzrank_taxes)) %>%
  print(include.rownames = TRUE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "kuz_rank_table.tex"))

#####--------------------------Table S39-----------------------------------#####
combined_trust %>%
  group_by(experiment) %>%
  summarise(pass_rate = mean(attention_heartattack, na.rm = TRUE))

table(combined_trust$attention_math)

# First stage analysis
study5_df <- 
  combined_trust %>% 
  filter(experiment == "MTurk_Pols17") %>%
  data.frame()

anova(lm(fmptrust_scale ~ Z + attention_math, data = study5_df),
      lm(fmptrust_scale ~ Z*attention_math, data = study5_df))

# ITT on federal spending items. 
anova(lm(fedspend_scale ~ Z + attention_math, data = study5_df),
      lm(fedspend_scale~ Z*attention_math, data = study5_df))


# 2SLS Estimates for both samples:
late_est_subset <- 
  iv_robust(as.formula(paste("fedspend_scale ~ fmptrust_scale + ", 
                             paste(X, collapse = "+"),
                             "| Z1 + ", paste(X, collapse = "+"))),
            data = study5_df  %>% filter(attention_math == 1)) %>% 
  tidy() %>%
  mutate_if(is.numeric, round, 2) %>% 
  filter(term == "fmptrust_scale") %>% 
  data.frame()

late_est_full <- 
  iv_robust(as.formula(paste("fedspend_scale ~ fmptrust_scale + ", 
                             paste(X, collapse = "+"),
                             "| Z1 + ", paste(X, collapse = "+"))),
            data = study5_df) %>% 
  tidy() %>%
  mutate_if(is.numeric, round, 2) %>% 
  filter(term == "fmptrust_scale") %>% 
  data.frame()

trust_est_subset <- 
  lm_robust(as.formula(paste("fmptrust_scale", "~ Z1 +", 
                             paste(covariates, collapse = "+"))),
            data = study5_df %>% filter(attention_math == 1)) %>%
  tidy() %>%
  mutate_if(is.numeric, round, 2) %>% 
  filter(term == "Z1") %>% 
  data.frame()

trust_est_full <- 
  lm_robust(as.formula(paste("fmptrust_scale", "~ Z1 +", 
                             paste(covariates, collapse = "+"))),
            data = study5_df) %>%
  tidy() %>%
  mutate_if(is.numeric, round, 2) %>% 
  filter(term == "Z1") %>% 
  data.frame()

redist_est_subset <- 
  lm_robust(as.formula(paste("fedspend_scale", "~ Z1 +", 
                             paste(covariates, collapse = "+"))),
            data = study5_df %>% filter(attention_math == 1)) %>%
  tidy() %>%
  mutate_if(is.numeric, round, 2) %>% 
  filter(term == "Z1") %>% 
  data.frame()

redist_est_full <- 
  lm_robust(as.formula(paste("fedspend_scale", "~ Z1 +", 
                             paste(covariates, collapse = "+"))),
            data = study5_df) %>%
  tidy() %>%
  mutate_if(is.numeric, round, 2) %>% 
  filter(term == "Z1") %>% 
  data.frame()


est_compare <- 
  bind_rows(trust_est_full, trust_est_subset, redist_est_full, 
            redist_est_subset, late_est_full, late_est_subset) %>% 
  mutate(outcome = ifelse(outcome == "fmptrust_scale", "Trust",
                          "Redistribution"),
         outcome = factor(outcome, levels = c("Trust", "Redistribution")),
         sample = rep(c(rep("Full", 1), rep("Subset", 1)), 3))
est_compare$term <- c("First Stage", "First Stage", "Reduced Form", 
                      "Reduced Form", "2SLS", "2SLS")

est_regtable <-
  est_compare %>%
  mutate(entry = make_entry(est = estimate, se = std.error, p = p.value)) %>% 
  dplyr::select(term, entry, outcome, sample) %>%
  spread(sample, entry) %>% 
  mutate(term = factor(term, 
                       levels = c("First Stage", "Reduced Form", "2SLS"))) %>% 
  arrange(term) %>% 
  select(-outcome)


xtable(est_regtable) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "robust_attention_check.tex"))

#####--------------------------Tables S40-S41------------------------------#####

# Op-Ed engagement/comprehension questions:
tmp <- 
  combined_trust %>%
  filter(experiment == "MTurk_Pols17") 

means <- 
  tmp %>%
  group_by(Z) %>%
  summarise(mean_est = mean(post_convince, na.rm = TRUE)) %>%
  mutate(measure = "Convincing") %>% 
  rbind(tmp %>%
          group_by(Z) %>% 
          summarise(mean_est = mean(post_access, na.rm = TRUE)) %>% 
          mutate(measure = "Accessible")) %>% 
  rbind(tmp %>%
          group_by(Z) %>% 
          summarise(mean_est = mean(post_priorinterest, na.rm = TRUE)) %>% 
          mutate(measure = "Prior Interest")) %>%
  rbind(tmp %>%
          group_by(Z) %>% 
          summarise(mean_est = mean(post_research, na.rm = TRUE)) %>% 
          mutate(measure = "Own Research"))

ses <- 
  tmp %>%
  filter(experiment == "MTurk_Pols17") %>%
  group_by(Z) %>%
  summarise(se = se_mean(post_convince, na.rm = TRUE)) %>%
  mutate(measure = "Convincing") %>% 
  rbind(combined_trust %>%
          group_by(Z) %>% 
          summarise(se = se_mean(post_access, na.rm = TRUE)) %>% 
          mutate(measure = "Accessible")) %>%
  rbind(tmp %>%
          group_by(Z) %>% 
          summarise(se = se_mean(post_priorinterest, na.rm = TRUE)) %>% 
          mutate(measure = "Prior Interest")) %>%
  rbind(tmp %>%
          group_by(Z) %>% 
          summarise(se = se_mean(post_research, na.rm = TRUE)) %>% 
          mutate(measure = "Own Research"))

n <- 
  combined_trust %>%
  group_by(Z) %>%
  summarise(n = sum(!is.na(post_convince)))

gg_df <- 
  left_join(n, left_join(means, ses)) %>%
  ungroup() %>%
  mutate(li = mean_est - 1.96*se,
         ui = mean_est + 1.96*se,
         Z = factor(Z, levels = c("Corrupt", "Placebo", "Honest"))) %>%
  data.frame

group_means <-
  gg_df %>%
  mutate(entry = table_entry(est = mean_est, se = se, digits = 2)) %>%
  select(measure, Z, entry) %>%
  spread(Z, entry)

xtable(group_means) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "post_oped_means.tex"))

# Visualization engagement/comprehension questions:
means <- 
  tmp %>%
  group_by(Z) %>%
  summarise(mean_est = mean(post_direction, na.rm = TRUE)) %>%
  mutate(measure = "Trend increased") %>% 
  rbind(tmp %>%
          group_by(Z) %>% 
          summarise(mean_est = mean(post_support, na.rm = TRUE)) %>% 
          mutate(measure = "Evidence supports")) 
ses <- 
  tmp %>%
  filter(experiment == "MTurk_Pols17") %>%
  group_by(Z) %>%
  summarise(se = se_mean(post_direction, na.rm = TRUE)) %>%
  mutate(measure = "Trend increased") %>% 
  rbind(combined_trust %>%
          group_by(Z) %>% 
          summarise(se = se_mean(post_support, na.rm = TRUE)) %>% 
          mutate(measure = "Evidence supports")) 
n <- 
  combined_trust %>%
  group_by(Z) %>%
  summarise(n = sum(!is.na(post_direction)))

gg_df <- 
  left_join(n, left_join(means, ses)) %>%
  ungroup() %>%
  mutate(li = mean_est - 1.96*se,
         ui = mean_est + 1.96*se,
         Z = factor(Z, levels = c("Corrupt", "Placebo", "Honest"))) %>%
  data.frame

group_means <-
  gg_df %>%
  mutate(entry = table_entry(est = mean_est, se = se, digits = 2)) %>%
  select(measure, Z, entry) %>%
  spread(Z, entry)

xtable(group_means) %>%
  print(include.rownames = FALSE, 
        include.colnames = FALSE,
        hline.after = c(),
        only.contents = TRUE,
        type = "latex",
        file = paste0(tabpath, "post_dataviz_means.tex"))
